#
# Replication file for
# "Identity, Industry, and Perceptions of Climate Futures," Journal of Politics
#
# Software: R 4.3.2 / macOS 14.2.1
# Processor: 2.4 GHz 8-Core Intel Core i9
# Memory: 64 GB
#
# Updated March 8, 2024
#

rm(list = ls())

#install.packages("pacman") # 0.5.1
pacman::p_load(tidyverse, # 2.0.0
               data.table, # 1.14.8
               DeclareDesign, # 1.0.4
               estimatr, # 1.0.0
               fixest, # 0.11.2
               ggpubr, # 0.6.0
               lfe, # 2.9-0
               magrittr, # 2.0.3
               mediation, # 4.5.0
               ri2, # 0.4.0
               stm, # 1.3.7
               tm) # 0.7-11

set.seed(1)

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# > Defining interaction plot function ####

interaction_plot_binary <- function(model, 
                                    effect,
                                    moderator, 
                                    interaction, 
                                    varcov="default", 
                                    conf=.95, 
                                    title="Marginal effects plot",
                                    xlabel="Value of moderator", 
                                    ylabel="Estimated marginal coefficient", 
                                    factor_labels=c(0,1)){
  
  # Extract Variance Covariance matrix
  if (varcov == "default"){
    covMat = vcov(model)
  }else{
    covMat = varcov
  }
  
  # Extract the data frame of the model
  mod_frame = model.frame(model)
  
  # Get coefficients of variables
  beta_1 = model$coefficients[[effect]]
  beta_3 = model$coefficients[[interaction]]
  
  # Create list of moderator values at which marginal effect is evaluated
  x_2 <- c(0,1)
  
  # Compute marginal effects
  delta_1 = beta_1 + beta_3*x_2
  
  # Compute variances
  var_1 = covMat[effect,effect] + (x_2^2)*covMat[interaction, interaction] + 2*x_2*covMat[effect, interaction]
  
  # Standard errors
  se_1 = sqrt(var_1)
  
  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score*se_1
  lower_bound = delta_1 - z_score*se_1
  
  # Determine the bounds of the graphing area
  max_y = max(upper_bound)
  min_y = min(lower_bound)
  
  # Initialize plotting window
  plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(-.5, 1.5), xlab=xlabel, ylab=ylabel, main=title, xaxt="n")
  
  # Plot points of estimated effects
  points(x=x_2, y=delta_1, pch=16)
  
  # Plot lines of confidence intervals
  lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1)
  points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg="black")
  lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1)
  points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg="black")
  
  # Label the axis
  axis(side=1, at=c(0,1), labels=factor_labels)
  
  # Add a dashed horizontal line for zero
  abline(h=0, lty=3)
  
}

# ****************** #
#################### #
#### DATA LOADING ####
#################### #
# ****************** #

# > Loading Lucid results ####

lucid <- read_rds("lucid1.rds")
lucid2 <- read_rds("lucid2.rds")

# > Loading QWI ####

qwi <- fread("qwi.csv")

# ************** #
################ #
#### ANALYSES ####
################ #
# ************** #

# > Figures 1-2. QWI descriptive plots ####

qwi21 <- qwi[year == 2019 & quarter == 4] # subsetting to 2019Q4

qwi21b <- qwi21[, .(WhiteEmp21 = sum(Emp[race == "A1" & ethnicity == "A1" & industry == "21"], na.rm = TRUE), # aggregating
                    BlackEmp21 = sum(Emp[race == "A2" & industry == "21"], na.rm = TRUE),
                    AllEmp21 = sum(Emp[industry == "21"], na.rm = TRUE),
                    WhiteEmp = sum(Emp[race == "A1" & ethnicity == "A1"], na.rm = TRUE),
                    BlackEmp = sum(Emp[race == "A2"], na.rm = TRUE),
                    AllEmp = sum(Emp, na.rm = TRUE)),
                by = geography]

qwi21b <- qwi21b[, WhiteEmpPct := WhiteEmp/AllEmp] # racial breakdown of county workforce
qwi21b <- qwi21b[, BlackEmpPct := BlackEmp/AllEmp]
qwi21b <- qwi21b[, NonWhiteEmpPct := 1 - WhiteEmp/AllEmp]

qwi21b <- qwi21b[, WhiteEmpPct21 := WhiteEmp21/AllEmp21] # racial breakdown of county mining workforce
qwi21b <- qwi21b[, BlackEmpPct21 := BlackEmp21/AllEmp21]
qwi21b <- qwi21b[, NonWhiteEmpPct21 := 1 - WhiteEmp21/AllEmp21]

qwi21b <- qwi21b[AllEmp21 >= 1] # limiting dataset to counties with at least one worker in mining

qwi21b <- qwi21b[, WhiteRel21 := WhiteEmpPct21 - WhiteEmpPct] # concentration in mining vs. overall
qwi21b <- qwi21b[, BlackRel21 := BlackEmpPct21 - BlackEmpPct]
qwi21b <- qwi21b[, NonWhiteRel21 := NonWhiteEmpPct21 - NonWhiteEmpPct]

# >> Climate-forcing ####

qwi21b %>%
  arrange(WhiteRel21) %>%
  mutate(id = 1:nrow(.)) %>%
  ggplot(aes(x = id)) +
  geom_hline(aes(yintercept = 0), lty = 2) +
  geom_point(aes(y = WhiteRel21, col = ifelse(WhiteRel21 > 0, "black", "gray"),
                 size = WhiteEmp21)) +
  scale_color_manual(values = c("black" = "black",
                                "gray" = "gray")) +
  theme_classic() +
  ylab("") +
  ggtitle("White Workforce Share") +
  theme(plot.title = element_text(size = 8, face = "bold"),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "none") +
  ylim(-1, 1) -> white_share_plot # Figure 1, left-hand panel

qwi21b %>%
  arrange(BlackRel21) %>%
  mutate(id = 1:nrow(.)) %>%
  ggplot(aes(x = id)) +
  geom_hline(aes(yintercept = 0), lty = 2) +
  geom_point(aes(y = BlackRel21, col = ifelse(BlackRel21 > 0, "black", "gray"),
                 size = BlackEmp21)) +
  scale_color_manual(values = c("black" = "black",
                                "gray" = "gray")) +
  theme_classic() +
  ylab("") +
  ggtitle("Black Workforce Share") +
  theme(plot.title = element_text(size = 8, face = "bold"),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "none") +
  ylim(-1, 1) -> black_share_plot  # Figure 1, center panel

qwi21b %>%
  arrange(NonWhiteRel21) %>%
  mutate(id = 1:nrow(.)) %>%
  ggplot(aes(x = id)) +
  geom_hline(aes(yintercept = 0), lty = 2) +
  geom_point(aes(y = NonWhiteRel21, col = ifelse(NonWhiteRel21 > 0, "black", "gray"),
                 size = (AllEmp21 - WhiteEmp21))) +
  scale_color_manual(values = c("black" = "black",
                                "gray" = "gray")) +
  theme_classic() +
  ylab("") +
  ggtitle("Minority Workforce Share") +
  theme(plot.title = element_text(size = 8, face = "bold"),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "none") +
  ylim(-1, 1) -> nonwhite_share_plot  # Figure 1, right-hand panel

ggarrange(white_share_plot,
          black_share_plot,
          nonwhite_share_plot,
          nrow = 1) # Figure 1 with all three panels

# >>> Appendix A, Figure A1 (by state) ####

qwi21b %>%
  group_by(substr(geography, 1, 2)) %>% # aggregating by state
  summarize(white = weighted.mean(WhiteRel21, w = WhiteEmp21),
            black = weighted.mean(BlackRel21, w = BlackEmp21),
            nonwhite = weighted.mean(NonWhiteRel21, w = (AllEmp21 - WhiteEmp21)),
            all_n = sum(AllEmp21, na.rm = TRUE),
            white_n = sum(WhiteEmp21, na.rm = TRUE),
            black_n = sum(BlackEmp21, na.rm = TRUE),
            nonwhite_n = sum((AllEmp21 - WhiteEmp21), na.rm = TRUE)) %>%
  ungroup() -> cf_skew_by_state

cf_skew_by_state %>%
  arrange(all_n) %>%
  mutate(id = 1:nrow(.)) %>%
  ggplot(aes(x = id)) +
  geom_hline(aes(yintercept = 0), lty = 2) +
  geom_point(aes(y = white,
                 size = white_n), col = "black") +
  scale_color_manual(values = c("black" = "black",
                                "red" = "red")) +
  theme_classic() +
  ylab("Skew Towards White Workers") +
  xlab("States in Ascending Order by Climate-Forcing Workforce Size (2019Q4)") +
  ggtitle("White Climate-Forcing Workforce Share by State") +
  theme(plot.title = element_text(size = 8, face = "bold"),
        axis.title.y = element_text(size = 8),
        axis.text.x = element_blank(),
        axis.title.x = element_text(size = 8),
        axis.ticks.x = element_blank(),
        legend.position = "none") +
  ylim(-.15, .35) -> white_share_plot_by_state # Appendix A, Figure A1, left-hand panel

# >> Climate-vulnerable ####

qwi11 <- qwi[year == 2019 & quarter == 4] # subsetting to 2019Q4

qwi11b <- qwi11[, .(WhiteEmp11 = sum(Emp[race == "A1" & ethnicity == "A1" & industry == "11"], na.rm = TRUE), # aggregating, in part filtered to agriculture, forestry, fishing, and hunting sector
                    BlackEmp11 = sum(Emp[race == "A2" & industry == "11"], na.rm = TRUE),
                    AllEmp11 = sum(Emp[industry == "11"], na.rm = TRUE),
                    WhiteEmp = sum(Emp[race == "A1" & ethnicity == "A1"], na.rm = TRUE),
                    BlackEmp = sum(Emp[race == "A2"], na.rm = TRUE),
                    AllEmp = sum(Emp, na.rm = TRUE)),
                by = geography]
qwi11b <- qwi11b[, WhiteEmpPct := WhiteEmp/AllEmp] # racial breakdowns of county workforce
qwi11b <- qwi11b[, BlackEmpPct := BlackEmp/AllEmp]
qwi11b <- qwi11b[, NonWhiteEmpPct := 1 - WhiteEmp/AllEmp]

qwi11b <- qwi11b[, WhiteEmpPct11 := WhiteEmp11/AllEmp11] # racial breakdowns of agriculture etc. sector
qwi11b <- qwi11b[, BlackEmpPct11 := BlackEmp11/AllEmp11]
qwi11b <- qwi11b[, NonWhiteEmpPct11 := 1 - WhiteEmp11/AllEmp11]

qwi11b <- qwi11b[AllEmp11 >= 1] # filtering to counties with at least one worker in agriculture etc.

qwi11b <- qwi11b[, WhiteRel11 := WhiteEmpPct11 - WhiteEmpPct]
qwi11b <- qwi11b[, BlackRel11 := BlackEmpPct11 - BlackEmpPct]
qwi11b <- qwi11b[, NonWhiteRel11 := NonWhiteEmpPct11 - NonWhiteEmpPct]

qwi11b %>%
  arrange(WhiteRel11) %>%
  mutate(id = 1:nrow(.)) %>%
  ggplot(aes(x = id)) +
  geom_hline(aes(yintercept = 0), lty = 2) +
  geom_point(aes(y = WhiteRel11, col = ifelse(WhiteRel11 > 0, "black", "gray"),
                 size = WhiteEmp11)) +
  scale_color_manual(values = c("black" = "black",
                                "gray" = "gray")) +
  theme_classic() +
  ylab("") +
  ggtitle("White Workforce Share") +
  theme(plot.title = element_text(size = 8, face = "bold"),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "none") +
  ylim(-1, 1) -> white_share_plot # Figure 2, left-hand panel

qwi11b %>%
  arrange(BlackRel11) %>%
  mutate(id = 1:nrow(.)) %>%
  ggplot(aes(x = id)) +
  geom_hline(aes(yintercept = 0), lty = 2) +
  geom_point(aes(y = BlackRel11, col = ifelse(BlackRel11 > 0, "black", "gray"),
                 size = BlackEmp11)) +
  scale_color_manual(values = c("black" = "black",
                                "gray" = "gray")) +
  theme_classic() +
  ylab("") +
  ggtitle("Black Workforce Share") +
  theme(plot.title = element_text(size = 8, face = "bold"),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "none") +
  ylim(-1, 1) -> black_share_plot # Figure 2, center panel

qwi11b %>%
  arrange(NonWhiteRel11) %>%
  mutate(id = 1:nrow(.)) %>%
  ggplot(aes(x = id)) +
  geom_hline(aes(yintercept = 0), lty = 2) +
  geom_point(aes(y = NonWhiteRel11, col = ifelse(NonWhiteRel11 > 0, "black", "gray"),
                 size = (AllEmp11 - WhiteEmp11))) +
  scale_color_manual(values = c("black" = "black",
                                "gray" = "gray")) +
  theme_classic() +
  ylab("") +
  ggtitle("Minority Workforce Share") +
  theme(plot.title = element_text(size = 8, face = "bold"),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "none") +
  ylim(-1, 1) -> nonwhite_share_plot # Figure 2, right-hand panel

ggarrange(white_share_plot,
          black_share_plot,
          nonwhite_share_plot,
          nrow = 1)  # Figure 2 with all three panels

# >>> Appendix A, Figure A1 (by state) ####

qwi11b %>%
  group_by(substr(geography, 1, 2)) %>%
  summarize(white = weighted.mean(WhiteRel11, w = WhiteEmp11),
            black = weighted.mean(BlackRel11, w = BlackEmp11),
            nonwhite = weighted.mean(NonWhiteRel11, w = (AllEmp11 - WhiteEmp11)),
            all_n = sum(AllEmp11, na.rm = TRUE),
            white_n = sum(WhiteEmp11, na.rm = TRUE),
            black_n = sum(BlackEmp11, na.rm = TRUE),
            nonwhite_n = sum((AllEmp11 - WhiteEmp11), na.rm = TRUE)) %>%
  ungroup() -> cv_skew_by_state

cv_skew_by_state %>%
  arrange(all_n) %>%
  mutate(id = 1:nrow(.)) %>%
  ggplot(aes(x = id)) +
  geom_hline(aes(yintercept = 0), lty = 2) +
  geom_point(aes(y = white,
                 size = white_n), col = "black") +
  scale_color_manual(values = c("black" = "black",
                                "red" = "red")) +
  theme_classic() +
  ylab("Skew Towards White Workers") +
  xlab("States in Ascending Order by Climate-Vulnerable Workforce Size (2019Q4)") +
  ggtitle("White Climate-Vulnerable Workforce Share by State") +
  theme(plot.title = element_text(size = 8, face = "bold"),
        axis.title.y = element_text(size = 8),
        axis.text.x = element_blank(),
        axis.title.x = element_text(size = 8),
        axis.ticks.x = element_blank(),
        legend.position = "none") +
  ylim(-.2, .1) -> white_share_plot_by_state_cv # Appendix A, Figure A1 (right-hand panel)

ggarrange(white_share_plot_by_state,
          white_share_plot_by_state_cv,
          nrow = 1) # Appendix A, Figure A1, with both panels

# > Figure 3. Control group beliefs ####

lucid %>%
  filter(treatment == "control") %>%
  summarize(decline_expectation_ff = mean(decline_expectation_ff2,
                                          na.rm = TRUE),
            decline_expectation_cv = mean(decline_expectation_cv2,
                                          na.rm = TRUE)) %>%
  ggplot() +
  xlim(-.3, .3) +
  geom_vline(aes(xintercept = 0), lty = 3) +
  geom_point(aes(y = "White-majority", x = decline_expectation_cv - decline_expectation_ff),
             color = "red") +
  geom_segment(aes(y = "White-majority", yend = "White-majority",
                   x = 0, xend = decline_expectation_cv - decline_expectation_ff),
               lty = 1,
               color = "red") +
  theme_minimal() +
  ggplot2::annotate("text", x = -.3, y = 1, label = "Climate-forcing\n industry", hjust = 0) +
  ggplot2::annotate("text", x = .3, y = 1, label = "Climate-vulnerable\n industry", hjust = 1) +
  scale_y_discrete(position = "bottom") +
  ggtitle("Perceived Risk of Decline") +
  ylab("") +
  xlab("Balance of Perceived Risk of Decline") -> decline_control_plot # Figure 3, left-hand panel
lucid %>%
  filter(treatment == "control") %$%
  t.test(decline_expectation_ff2, decline_expectation_cv2) # t-test comparison of decline expectations for climate-forcing vs. climate-vulnerable industries

lucid %>%
  filter(treatment == "control") %>%
  summarize(subsidy_ff = mean(subsidy_expectation_ff2,
                              na.rm = TRUE)) %>%
  ggplot() +
  xlim(-.3, .3) +
  geom_vline(aes(xintercept = 0), lty = 3) +
  geom_point(aes(y = "White-majority", x = .5 - subsidy_ff),
             color = "blue") +
  geom_segment(aes(y = "White-majority", yend = "White-majority",
                   x = 0, xend = .5 - subsidy_ff),
               lty = 1,
               color = "blue") +
  theme_minimal() +
  ggplot2::annotate("text", x = -.3, y = 1, label = "Climate-forcing\nindustry", hjust = 0) +
  ggplot2::annotate("text", x = .3, y = 1, label = "Climate-vulnerable\nindustry", hjust = 1) +
  scale_y_discrete(position = "bottom") +
  ggtitle("Subsidy Expectations") +
  ylab("") +
  xlab("Balance of Perceived Likelihood of Subsidization") -> subsidy_control_plot # Figure 3, right-hand panel

ggarrange(decline_control_plot,
          subsidy_control_plot,
          nrow = 1) # Figure 3

# > Perceptions of climate risk ####

# >> Figure 4. By prior perception of bias ####

lucid %>%
  filter(treatment != "control") %>%
  lm_robust(decline_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg1a # Figure 4, Model 1
lucid %>%
  filter(treatment != "control") %>%
  lm_robust(decline_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg1  # Figure 4, Model 2

lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(decline_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg2a # Figure 4, Model 3
lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(decline_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg2 # Figure 4, Model 4

lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 0) %>%
  lm_robust(decline_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg3a # Figure 4, Model 5
lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 0) %>%
  lm_robust(decline_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg3  # Figure 4, Model 6

#NOTE: the coefficients in Figure 4 are inverse of what's shown here for Models 7-12,
#as "treatment" indicates a Black-majority climate-vulnerable industry,
#while the table lists results for the inverse (white-majority climate-vulnerable industry)
lucid %>%
  filter(treatment != "control") %>%
  lm_robust(decline_expectation_cv ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg4a # Figure 4, Model 7
lucid %>%
  filter(treatment != "control") %>%
  lm_robust(decline_expectation_cv ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg4 # Figure 4, Model 8

lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(decline_expectation_cv ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg5a # Figure 4, Model 9
lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(decline_expectation_cv ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg5 # Figure 4, Model 10

lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 0) %>%
  lm_robust(decline_expectation_cv ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg6a # Figure 4, Model 11
lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 0) %>%
  lm_robust(decline_expectation_cv ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg6 # Figure 4, Model 12

lucid %>%
  filter(favW_tok == 0) %>%
  summarize(decline_ff_w = mean(decline_expectation_ff[treatment == "ffWhite_cvBlack"],
                                na.rm = TRUE),
            decline_cv_b = mean(decline_expectation_cv[treatment == "ffWhite_cvBlack"],
                                na.rm = TRUE),
            decline_ff_b = mean(decline_expectation_ff[treatment == "ffBlack_cvWhite"],
                                na.rm = TRUE),
            decline_cv_w = mean(decline_expectation_cv[treatment == "ffBlack_cvWhite"],
                                na.rm = TRUE),
            decline_ff_c = mean(decline_expectation_ff2[treatment == "control"],
                                na.rm = TRUE),
            decline_cv_c = mean(decline_expectation_cv2[treatment == "control"],
                                na.rm = TRUE)) %>%
  ggplot() +
  xlim(-1.2, 1.2) +
  geom_vline(aes(xintercept = 0), lty = 3) +
  geom_point(aes(y = "White-majority", x = decline_cv_b - decline_ff_w),
             color = "red") +
  geom_segment(aes(y = "White-majority", yend = "White-majority",
                   x = 0, xend = decline_cv_b - decline_ff_w),
               lty = 1,
               color = "red") +
  geom_point(aes(y = "Control", x = decline_cv_c - decline_ff_c),
             color = "gray50") +
  geom_segment(aes(y = "Control", yend = "Control", 
                   x = 0, xend = decline_cv_c - decline_ff_c),
               lty = 1,
               color = "gray50") +
  geom_point(aes(y = "Black-majority", x = decline_cv_w - decline_ff_b),
             color = "red") +
  geom_segment(aes(y = "Black-majority", yend = "Black-majority", 
                   x = 0, xend = decline_cv_w - decline_ff_b),
               lty = 1,
               color = "red") +
  theme_minimal() +
  ggplot2::annotate("text", x = -1.2, y = 3, label = "White-maj.\nCF industry", hjust = 0) +
  ggplot2::annotate("text", x = 1.2, y = 3, label = "Black-maj.\nCV industry", hjust = 1) +
  ggplot2::annotate("text", x = -1.2, y = 2, label = "CF industry\n(control)", hjust = 0) +
  ggplot2::annotate("text", x = 1.2, y = 2, label = "CV industry\n(control)", hjust = 1) +
  ggplot2::annotate("text", x = -1.2, y = 1, label = "Black-maj.\nCF industry", hjust = 0) +
  ggplot2::annotate("text", x = 1.2, y = 1, label = "White-maj.\nCV industry", hjust = 1) +
  scale_y_discrete(position = "bottom") +
  ggtitle("State Does Not Favor White Citizens") +
  ylab("") +
  xlab("Balance of Perceived Risk of Decline") -> decline_not_fav_plot # Figure 4, right-hand graphic
lucid %>%
  filter(favW_tok == 1) %>%
  summarize(decline_ff_w = mean(decline_expectation_ff[treatment == "ffWhite_cvBlack"],
                                na.rm = TRUE),
            decline_cv_b = mean(decline_expectation_cv[treatment == "ffWhite_cvBlack"],
                                na.rm = TRUE),
            decline_ff_b = mean(decline_expectation_ff[treatment == "ffBlack_cvWhite"],
                                na.rm = TRUE),
            decline_cv_w = mean(decline_expectation_cv[treatment == "ffBlack_cvWhite"],
                                na.rm = TRUE),
            decline_ff_c = mean(decline_expectation_ff2[treatment == "control"],
                                na.rm = TRUE),
            decline_cv_c = mean(decline_expectation_cv2[treatment == "control"],
                                na.rm = TRUE)) %>%
  ggplot() +
  xlim(-1.2, 1.2) +
  geom_vline(aes(xintercept = 0), lty = 3) +
  geom_point(aes(y = "White-majority", x = decline_cv_b - decline_ff_w),
             color = "red") +
  geom_segment(aes(y = "White-majority", yend = "White-majority",
                   x = 0, xend = decline_cv_b - decline_ff_w),
               lty = 1,
               color = "red") +
  geom_point(aes(y = "Control", x = decline_cv_c - decline_ff_c),
             color = "gray50") +
  geom_segment(aes(y = "Control", yend = "Control", 
                   x = 0, xend = decline_cv_c - decline_ff_c),
               lty = 1,
               color = "gray50") +
  geom_point(aes(y = "Black-majority", x = decline_cv_w - decline_ff_b),
             color = "red") +
  geom_segment(aes(y = "Black-majority", yend = "Black-majority", 
                   x = 0, xend = decline_cv_w - decline_ff_b),
               lty = 1,
               color = "red") +
  theme_minimal() +
  ggplot2::annotate("text", x = -1.2, y = 3, label = "White-maj.\nCF industry", hjust = 0) +
  ggplot2::annotate("text", x = 1.2, y = 3, label = "Black-maj.\nCV industry", hjust = 1) +
  ggplot2::annotate("text", x = -1.2, y = 2, label = "CF industry\n(control)", hjust = 0) +
  ggplot2::annotate("text", x = 1.2, y = 2, label = "CV industry\n(control)", hjust = 1) +
  ggplot2::annotate("text", x = -1.2, y = 1, label = "Black-maj.\nCF industry", hjust = 0) +
  ggplot2::annotate("text", x = 1.2, y = 1, label = "White-maj.\nCV industry", hjust = 1) +
  scale_y_discrete(position = "bottom") +
  ggtitle("State Favors White Citizens") +
  ylab("") +
  xlab("Balance of Perceived Risk of Decline") -> decline_fav_plot  # Figure 4, left-hand graphic

ggarrange(decline_fav_plot,
          decline_not_fav_plot,
          nrow = 1) # Figure 4 graphic

# >> Footnote 30 ####
lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(I(decline_expectation_ff - decline_expectation_cv) ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE)

# > Access to subsidies ####

# >> Figure 5. By prior perception of bias ####

lucid %>%
  filter(treatment != "control") %>%
  lm_robust(subsidy_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg1a # Figure 5, Model 1
lucid %>%
  filter(treatment != "control") %>%
  lm_robust(subsidy_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg1 # Figure 5, Model 2

lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(subsidy_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg2a  # Figure 5, Model 3
lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(subsidy_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg2  # Figure 5, Model 4

lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 0) %>%
  lm_robust(subsidy_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg3a  # Figure 5, Model 5
lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 0) %>%
  lm_robust(subsidy_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg3  # Figure 5, Model 6

lucid %>%
  filter(favW_tok == 0) %>%
  summarize(subsidy_ff_w = mean(subsidy_expectation_ff[treatment == "ffWhite_cvBlack"],
                                na.rm = TRUE),
            subsidy_ff_b = mean(subsidy_expectation_ff[treatment == "ffBlack_cvWhite"],
                                na.rm = TRUE),
            subsidy_ff_c = mean(subsidy_expectation_ff2[treatment == "control"],
                                na.rm = TRUE)) %>%
  ggplot() +
  xlim(-.3, .3) +
  geom_vline(aes(xintercept = 0), lty = 3) +
  geom_point(aes(y = "White-majority", x = .5 - subsidy_ff_w),
             color = "blue") +
  geom_segment(aes(y = "White-majority", yend = "White-majority",
                   x = 0, xend = .5 - subsidy_ff_w),
               lty = 2,
               color = "blue") +
  geom_point(aes(y = "Control", x = .5 - subsidy_ff_c),
             color = "gray50") +
  geom_segment(aes(y = "Control", yend = "Control",
                   x = 0, xend = .5 - subsidy_ff_c),
               lty = 2,
               color = "gray50") +
  geom_point(aes(y = "Black-majority", x = .5 - subsidy_ff_b),
             color = "blue") +
  geom_segment(aes(y = "Black-majority", yend = "Black-majority", 
                   x = 0, xend = .5 - subsidy_ff_b),
               lty = 2,
               color = "blue") +
  theme_minimal() +
  ggplot2::annotate("text", x = -.3, y = 3, label = "White-maj.\nCF industry", hjust = 0) +
  ggplot2::annotate("text", x = .3, y = 3, label = "Black-maj.\nCV industry", hjust = 1) +
  ggplot2::annotate("text", x = -.3, y = 2, label = "CF industry\n(control)", hjust = 0) +
  ggplot2::annotate("text", x = .3, y = 2, label = "CV industry\n(control)", hjust = 1) +
  ggplot2::annotate("text", x = -.3, y = 1, label = "Black-maj.\nCF industry", hjust = 0) +
  ggplot2::annotate("text", x = .3, y = 1, label = "White-maj.\nCV industry", hjust = 1) +
  scale_y_discrete(position = "bottom") +
  ggtitle("State Does Not Favor White Citizens") +
  ylab("") +
  xlab("Balance of Perceived Likelihood of Subsidization") -> subsidy_not_fav_plot # Figure 5 right-hand graphic
lucid %>%
  filter(favW_tok == 1) %>%
  summarize(subsidy_ff_w = mean(subsidy_expectation_ff[treatment == "ffWhite_cvBlack"],
                                na.rm = TRUE),
            subsidy_ff_b = mean(subsidy_expectation_ff[treatment == "ffBlack_cvWhite"],
                                na.rm = TRUE),
            subsidy_ff_c = mean(subsidy_expectation_ff2[treatment == "control"],
                                na.rm = TRUE)) %>%
  ggplot() +
  xlim(-.3, .3) +
  geom_vline(aes(xintercept = 0), lty = 3) +
  geom_point(aes(y = "White-majority", x = .5 - subsidy_ff_w),
             color = "blue") +
  geom_segment(aes(y = "White-majority", yend = "White-majority",
                   x = 0, xend = .5 - subsidy_ff_w),
               lty = 2,
               color = "blue") +
  geom_point(aes(y = "Control", x = .5 - subsidy_ff_c),
             color = "gray50") +
  geom_segment(aes(y = "Control", yend = "Control",
                   x = 0, xend = .5 - subsidy_ff_c),
               lty = 2,
               color = "gray50") +
  geom_point(aes(y = "Black-majority", x = .5 - subsidy_ff_b),
             color = "blue") +
  geom_segment(aes(y = "Black-majority", yend = "Black-majority", 
                   x = 0, xend = .5 - subsidy_ff_b),
               lty = 2,
               color = "blue") +
  theme_minimal() +
  ggplot2::annotate("text", x = -.3, y = 3, label = "White-maj.\nCF industry", hjust = 0) +
  ggplot2::annotate("text", x = .3, y = 3, label = "Black-maj.\nCV industry", hjust = 1) +
  ggplot2::annotate("text", x = -.3, y = 2, label = "CF industry\n(control)", hjust = 0) +
  ggplot2::annotate("text", x = .3, y = 2, label = "CV industry\n(control)", hjust = 1) +
  ggplot2::annotate("text", x = -.3, y = 1, label = "Black-maj.\nCF industry", hjust = 0) +
  ggplot2::annotate("text", x = .3, y = 1, label = "White-maj.\nCV industry", hjust = 1) +
  scale_y_discrete(position = "bottom") +
  ggtitle("State Favors White Citizens") +
  ylab("") +
  xlab("Balance of Perceived Likelihood of Subsidization") -> subsidy_fav_plot  # Figure 5 left-hand graphic

ggarrange(subsidy_fav_plot,
          subsidy_not_fav_plot,
          nrow = 1) # Figure 5 graphic

# >> Footnote 31. Decline expectations ~ subsidy expectations ####

lucid %>% # climate-forcing regression
  filter(treatment != "control") %>%
  lm_robust(decline_expectation_ff ~ subsidy_expectation_ff,
            data = .,
            fixed_effects = ~ NONHWHITE) %>%
  tidy()

lucid %>% # climate-vulnerable regression
  filter(treatment != "control") %>%
  lm_robust(decline_expectation_cv ~ subsidy_expectation_ff,
            data = .,
            fixed_effects = ~ NONHWHITE) %>%
  tidy()

# > Worker protection and political efficacy #### 

# >> Table 1. By prior perception of bias ####

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     na.action = na.omit) -> bias_protection_reg1 # Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1,
     na.action = na.omit) -> bias_protection_reg2 # Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0,
     na.action = na.omit) -> bias_protection_reg3 # Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     na.action = na.omit) -> bias_protection_reg4 # Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1,
     na.action = na.omit) -> bias_protection_reg5 # Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0,
     na.action = na.omit) -> bias_protection_reg6 # Model 6

# **************** #
################## #
#### APPENDICES ####
################## #
# **************** #

# > Appendix B. Balance ####

lucid %>%
  filter(!is.na(treatment)) %>%
  group_by(treatment) %>%
  summarize(female = mean(gender == "Female", na.rm = TRUE),
            age = mean(age, na.rm = TRUE),
            income = mean(income, na.rm = TRUE),
            nonhispwhite = mean(NONHWHITE, na.rm = TRUE),
            republican = mean(party_fct == "Republican", na.rm = TRUE),
            ideology = mean(ideology, na.rm = TRUE),
            education = mean(education, na.rm = TRUE)) %>%
  ungroup()

t.test(lucid$gender[lucid$treatment == "ffBlack_cvWhite"] == "Female",
       lucid$gender[lucid$treatment == "ffWhite_cvBlack"] == "Female")
t.test(lucid$age[lucid$treatment == "ffBlack_cvWhite"],
       lucid$age[lucid$treatment == "ffWhite_cvBlack"])
t.test(lucid$income[lucid$treatment == "ffBlack_cvWhite"],
       lucid$income[lucid$treatment == "ffWhite_cvBlack"])
t.test(lucid$NONHWHITE[lucid$treatment == "ffBlack_cvWhite"],
       lucid$NONHWHITE[lucid$treatment == "ffWhite_cvBlack"])
t.test(lucid$party_fct[lucid$treatment == "ffBlack_cvWhite"] == "Republican",
       lucid$party_fct[lucid$treatment == "ffWhite_cvBlack"] == "Republican")
t.test(lucid$ideology[lucid$treatment == "ffBlack_cvWhite"],
       lucid$ideology[lucid$treatment == "ffWhite_cvBlack"])
t.test(lucid$education[lucid$treatment == "ffBlack_cvWhite"],
       lucid$education[lucid$treatment == "ffWhite_cvBlack"])

# > Appendix D: Outcome Summary Statistics ####

lucid %>%
  filter(treatment == "control") %>%
  summarize_at(vars(decline_expectation_ff2), list(length, mean, sd, min, max)) # summary statistics for perceived risk of decline in climate-forcing industry
lucid %>%
  filter(treatment == "control") %>%
  filter(!is.na(decline_expectation_cv2)) %>%
  summarize_at(vars(decline_expectation_cv2), list(length, mean, sd, min, max)) # summary statistics for perceived risk of decline in climate-vulnerable industry
lucid %>%
  filter(treatment == "control") %>%
  filter(!is.na(subsidy_expectation_ff2)) %>%
  summarize_at(vars(subsidy_expectation_ff2), list(length, mean, sd, min, max)) # summary statistics for perceived probability of climate-forcing industry winning subsidies

# > Appendix E. Perceptions of Climate Risk: Bonferroni Correction ####

c(bias_reg1$p.value[[1]],
  bias_reg2$p.value[[1]],
  bias_reg3$p.value[[1]],
  bias_reg4$p.value[[1]],
  bias_reg5$p.value[[1]],
  bias_reg6$p.value[[1]]) -> risk_p_vals
p.adjust(p = risk_p_vals, method = "bonferroni")

# > Appendix F. Perceptions of Climate Risk: Randomization Inference ####

# >> All - FF ####

lucid_ri_fav <- lucid %>%
  dplyr::select(treatment, 
                decline_expectation_ff, 
                decline_expectation_cv,
                NONHWHITE,
                gender_fct,
                ideology,
                party_fct,
                age,
                income,
                education) %>%
  mutate(treatment = case_when(treatment == "ffWhite_cvBlack" ~ 1,
                               treatment == "ffBlack_cvWhite" ~ 0,
                               treatment == "control" ~ NA_real_,
                               TRUE ~ NA_real_)) %>%
  rename(Z = treatment,
         Yf = decline_expectation_ff,
         Yv = decline_expectation_cv,
         block = NONHWHITE) %>%
  filter(!is.na(Z) & !is.na(block) & !is.na(Yf) & !is.na(Yv) &
           !is.na(gender_fct) & !is.na(ideology) & !is.na(party_fct) &
           !is.na(age) & !is.na(income) & !is.na(education))

with(lucid_ri_fav, table(block, Z))

block_m <- with(lucid_ri_fav, tapply(Z, block, sum) / 2)

declaration <- with(lucid_ri_fav, {
  declare_ra(blocks = block,
             block_m = block_m)
})

riskF_ri_out <- conduct_ri(Yf ~ Z +
                             gender_fct + ideology + party_fct + age + income + education,
                           sharp_hypothesis = 0,
                           declaration = declaration,
                           data = lucid_ri_fav[!is.na(lucid_ri_fav$Z) &
                                                 !is.na(lucid_ri_fav$Yf), ])
tidy(riskF_ri_out) # Table F1, Model 1

# >> Prior: whites favored - FF ####

lucid_ri_fav2 <- lucid %>%
  filter(favW_tok == 1) %>%
  dplyr::select(treatment, 
                decline_expectation_ff, 
                decline_expectation_cv,
                NONHWHITE,
                gender_fct,
                ideology,
                party_fct,
                age,
                income,
                education) %>%
  mutate(treatment = case_when(treatment == "ffWhite_cvBlack" ~ 1,
                               treatment == "ffBlack_cvWhite" ~ 0,
                               treatment == "control" ~ NA_real_,
                               TRUE ~ NA_real_)) %>%
  rename(Z = treatment,
         Yf = decline_expectation_ff,
         Yv = decline_expectation_cv,
         block = NONHWHITE) %>%
  filter(!is.na(Z) & !is.na(block) & !is.na(Yf) & !is.na(Yv) &
           !is.na(gender_fct) & !is.na(ideology) & !is.na(party_fct) &
           !is.na(age) & !is.na(income) & !is.na(education))

with(lucid_ri_fav2, table(block, Z))

block_m2 <- with(lucid_ri_fav2, tapply(Z, block, sum) / 2)

declaration2 <- with(lucid_ri_fav2, {
  declare_ra(blocks = block,
             block_m = block_m2)
})

riskF_ri_out2 <- conduct_ri(Yf ~ Z +
                              gender_fct + ideology + party_fct + age + income + education,
                            sharp_hypothesis = 0,
                            declaration = declaration2,
                            data = lucid_ri_fav2[!is.na(lucid_ri_fav2$Z) &
                                                   !is.na(lucid_ri_fav2$Yf), ])
tidy(riskF_ri_out2) # Table F1, Model 2

# >> Prior: whites favored - FF ####

lucid_ri_fav3 <- lucid %>%
  filter(favW_tok == 0) %>%
  dplyr::select(treatment, 
                decline_expectation_ff, 
                decline_expectation_cv,
                NONHWHITE,
                gender_fct,
                ideology,
                party_fct,
                age,
                income,
                education) %>%
  mutate(treatment = case_when(treatment == "ffWhite_cvBlack" ~ 1,
                               treatment == "ffBlack_cvWhite" ~ 0,
                               treatment == "control" ~ NA_real_,
                               TRUE ~ NA_real_)) %>%
  rename(Z = treatment,
         Yf = decline_expectation_ff,
         Yv = decline_expectation_cv,
         block = NONHWHITE) %>%
  filter(!is.na(Z) & !is.na(block) & !is.na(Yf) & !is.na(Yv) &
           !is.na(gender_fct) & !is.na(ideology) & !is.na(party_fct) &
           !is.na(age) & !is.na(income) & !is.na(education))

with(lucid_ri_fav3, table(block, Z))

block_m3 <- with(lucid_ri_fav3, tapply(Z, block, sum) / 2)

declaration3 <- with(lucid_ri_fav3, {
  declare_ra(blocks = block,
             block_m = block_m3)
})

riskF_ri_out3 <- conduct_ri(Yf ~ Z +
                              gender_fct + ideology + party_fct + age + income + education,
                            sharp_hypothesis = 0,
                            declaration = declaration3,
                            data = lucid_ri_fav3[!is.na(lucid_ri_fav3$Z) &
                                                   !is.na(lucid_ri_fav3$Yf), ])
tidy(riskF_ri_out3) # Table F1, Model 3

# >> All - CV ####

lucid_ri_fav <- lucid %>%
  dplyr::select(treatment, 
                decline_expectation_ff, 
                decline_expectation_cv,
                NONHWHITE,
                gender_fct,
                ideology,
                party_fct,
                age,
                income,
                education) %>%
  mutate(treatment = case_when(treatment == "ffWhite_cvBlack" ~ 1,
                               treatment == "ffBlack_cvWhite" ~ 0,
                               treatment == "control" ~ NA_real_,
                               TRUE ~ NA_real_)) %>%
  rename(Z = treatment,
         Yf = decline_expectation_ff,
         Yv = decline_expectation_cv,
         block = NONHWHITE) %>%
  filter(!is.na(Z) & !is.na(block) & !is.na(Yf) & !is.na(Yv) &
           !is.na(gender_fct) & !is.na(ideology) & !is.na(party_fct) &
           !is.na(age) & !is.na(income) & !is.na(education))

with(lucid_ri_fav, table(block, Z))

block_m <- with(lucid_ri_fav, tapply(Z, block, sum) / 2)

declaration <- with(lucid_ri_fav, {
  declare_ra(blocks = block,
             block_m = block_m)
})

riskV_ri_out <- conduct_ri(Yv ~ Z +
                             gender_fct + ideology + party_fct + age + income + education,
                           sharp_hypothesis = 0,
                           declaration = declaration,
                           data = lucid_ri_fav[!is.na(lucid_ri_fav$Z) &
                                                 !is.na(lucid_ri_fav$Yv), ])
tidy(riskV_ri_out) # Table F1, Model 4

# >> Prior: whites favored - CV ####

lucid_ri_fav2 <- lucid %>%
  filter(favW_tok == 1) %>%
  dplyr::select(treatment, 
                decline_expectation_ff, 
                decline_expectation_cv,
                NONHWHITE,
                gender_fct,
                ideology,
                party_fct,
                age,
                income,
                education) %>%
  mutate(treatment = case_when(treatment == "ffWhite_cvBlack" ~ 1,
                               treatment == "ffBlack_cvWhite" ~ 0,
                               treatment == "control" ~ NA_real_,
                               TRUE ~ NA_real_)) %>%
  rename(Z = treatment,
         Yf = decline_expectation_ff,
         Yv = decline_expectation_cv,
         block = NONHWHITE) %>%
  filter(!is.na(Z) & !is.na(block) & !is.na(Yf) & !is.na(Yv) &
           !is.na(gender_fct) & !is.na(ideology) & !is.na(party_fct) &
           !is.na(age) & !is.na(income) & !is.na(education))

with(lucid_ri_fav2, table(block, Z))

block_m2 <- with(lucid_ri_fav2, tapply(Z, block, sum) / 2)

declaration2 <- with(lucid_ri_fav2, {
  declare_ra(blocks = block,
             block_m = block_m2)
})

riskV_ri_out2 <- conduct_ri(Yv ~ Z +
                              gender_fct + ideology + party_fct + age + income + education,
                            sharp_hypothesis = 0,
                            declaration = declaration2,
                            data = lucid_ri_fav2[!is.na(lucid_ri_fav2$Z) &
                                                   !is.na(lucid_ri_fav2$Yv), ])
tidy(riskV_ri_out2) # Table F1, Model 5

# >> Prior: whites favored - CV ####

lucid_ri_fav3 <- lucid %>%
  filter(favW_tok == 0) %>%
  dplyr::select(treatment, 
                decline_expectation_ff, 
                decline_expectation_cv,
                NONHWHITE,
                gender_fct,
                ideology,
                party_fct,
                age,
                income,
                education) %>%
  mutate(treatment = case_when(treatment == "ffWhite_cvBlack" ~ 1,
                               treatment == "ffBlack_cvWhite" ~ 0,
                               treatment == "control" ~ NA_real_,
                               TRUE ~ NA_real_)) %>%
  rename(Z = treatment,
         Yf = decline_expectation_ff,
         Yv = decline_expectation_cv,
         block = NONHWHITE) %>%
  filter(!is.na(Z) & !is.na(block) & !is.na(Yf) & !is.na(Yv) &
           !is.na(gender_fct) & !is.na(ideology) & !is.na(party_fct) &
           !is.na(age) & !is.na(income) & !is.na(education))

with(lucid_ri_fav3, table(block, Z))

block_m3 <- with(lucid_ri_fav3, tapply(Z, block, sum) / 2)

declaration3 <- with(lucid_ri_fav3, {
  declare_ra(blocks = block,
             block_m = block_m3)
})

riskV_ri_out3 <- conduct_ri(Yv ~ Z +
                              gender_fct + ideology + party_fct + age + income + education,
                            sharp_hypothesis = 0,
                            declaration = declaration3,
                            data = lucid_ri_fav3[!is.na(lucid_ri_fav3$Z) &
                                                   !is.na(lucid_ri_fav3$Yv), ])
tidy(riskV_ri_out3) # Table F1, Model 6

# > Appendix G. Perceptions of Climate Risk: Alternative Bias Measures ####

# >> G1. White people have an easier time... ####

lucid %>%
  filter(treatment != "control") %>%
  filter(whitesupport %in% 3:4) %>%
  lm_robust(decline_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg2a_alt1 # G1, Model 1
lucid %>%
  filter(treatment != "control") %>%
  filter(whitesupport %in% 3:4) %>%
  lm_robust(decline_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg2_alt1 # G1, Model 2

lucid %>%
  filter(treatment != "control") %>%
  filter(whitesupport %in% 0:2) %>%
  lm_robust(decline_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg3a_alt1 # G1, Model 3
lucid %>%
  filter(treatment != "control") %>%
  filter(whitesupport %in% 0:2) %>%
  lm_robust(decline_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg3_alt1 # G1, Model 4

#NOTE: the coefficients in G1 are inverse of what's shown here for Models 5-8,
#as "treatment" indicates a Black-majority climate-vulnerable industry,
#while the table lists results for the inverse (white-majority climate-vulnerable industry)
lucid %>%
  filter(treatment != "control") %>%
  filter(whitesupport %in% 3:4) %>%
  lm_robust(decline_expectation_cv ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg5a_alt1 # G1, Model 5
lucid %>%
  filter(treatment != "control") %>%
  filter(whitesupport %in% 3:4) %>%
  lm_robust(decline_expectation_cv ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg5_alt1 # G1, Model 6

lucid %>%
  filter(treatment != "control") %>%
  filter(whitesupport %in% 0:2) %>%
  lm_robust(decline_expectation_cv ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg6a_alt1 # G1, Model 7
lucid %>%
  filter(treatment != "control") %>%
  filter(whitesupport %in% 0:2) %>%
  felm(decline_expectation_cv ~ treatment  +
         gender_fct + ideology + party_fct + age + income + education | NONHWHITE | 0 | 0,
       data = .) -> bias_reg6_alt1 # G1, Model 8

# >> G2. On average, politicians are today more responsive... ####

lucid %>%
  filter(treatment != "control") %>%
  filter(whiteresponsiveness %in% 3:4) %>%
  lm_robust(decline_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg2a_alt2 # G2, Model 1
lucid %>%
  filter(treatment != "control") %>%
  filter(whiteresponsiveness %in% 3:4) %>%
  lm_robust(decline_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg2_alt2 # G2, Model 2

lucid %>%
  filter(treatment != "control") %>%
  filter(whiteresponsiveness %in% 0:2) %>%
  lm_robust(decline_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg3a_alt2 # G2, Model 3
lucid %>%
  filter(treatment != "control") %>%
  filter(whiteresponsiveness %in% 0:2) %>%
  lm_robust(decline_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg3_alt2 # G2, Model 4

#NOTE: the coefficients in G2 are inverse of what's shown here for Models 5-8,
#as "treatment" indicates a Black-majority climate-vulnerable industry,
#while the table lists results for the inverse (white-majority climate-vulnerable industry)
lucid %>%
  filter(treatment != "control") %>%
  filter(whiteresponsiveness %in% 3:4) %>%
  lm_robust(decline_expectation_cv ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg5a_alt2 # G2, Model 5
lucid %>%
  filter(treatment != "control") %>%
  filter(whiteresponsiveness %in% 3:4) %>%
  lm_robust(decline_expectation_cv ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg5_alt2 # G2, Model 6

lucid %>%
  filter(treatment != "control") %>%
  filter(whiteresponsiveness %in% 0:2) %>%
  lm_robust(decline_expectation_cv ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg6a_alt2 # G2, Model 7
lucid %>%
  filter(treatment != "control") %>%
  filter(whiteresponsiveness %in% 0:2) %>%
  lm_robust(decline_expectation_cv ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg6_alt2 # G2, Model 8

# >> G3. When the government spends money, white people... ####

lucid %>%
  filter(treatment != "control") %>%
  filter(whitemoneywins %in% 3:4) %>%
  lm_robust(decline_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg2a_alt3 # G3, Model 1
lucid %>%
  filter(treatment != "control") %>%
  filter(whitemoneywins %in% 3:4) %>%
  lm_robust(decline_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg2_alt3 # G3, Model 2

lucid %>%
  filter(treatment != "control") %>%
  filter(whitemoneywins %in% 0:2) %>%
  lm_robust(decline_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg3a_alt3 # G3, Model 3
lucid %>%
  filter(treatment != "control") %>%
  filter(whitemoneywins %in% 0:2) %>%
  lm_robust(decline_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg3_alt3 # G3, Model 4

#NOTE: the coefficients in G3 are inverse of what's shown here for Models 5-8,
#as "treatment" indicates a Black-majority climate-vulnerable industry,
#while the table lists results for the inverse (white-majority climate-vulnerable industry)
lucid %>%
  filter(treatment != "control") %>%
  filter(whitemoneywins %in% 3:4) %>%
  lm_robust(decline_expectation_cv ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg5a_alt3 # G3, Model 5
lucid %>%
  filter(treatment != "control") %>%
  filter(whitemoneywins %in% 3:4) %>%
  lm_robust(decline_expectation_cv ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg5_alt3 # G3, Model 6

lucid %>%
  filter(treatment != "control") %>%
  filter(whitemoneywins %in% 0:2) %>%
  lm_robust(decline_expectation_cv ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg6a_alt3 # G3, Model 7
lucid %>%
  filter(treatment != "control") %>%
  filter(whitemoneywins %in% 0:2) %>%
  lm_robust(decline_expectation_cv ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg6_alt3 # G3, Model 8

# > Appendix H. Perceptions of Climate Risk: Interaction Models ####

lucid %>%
  filter(treatment != "control") %>%
  lm(decline_expectation_ff ~ treatment*favW_tok +
       gender_fct + ideology + party_fct + age + income + education +
       as.factor(NONHWHITE),
     data = .) -> risk_intF
interaction_plot_binary(model = risk_intF,
                        effect = "treatmentffWhite_cvBlack",
                        moderator = "favW_tok", 
                        interaction = "treatmentffWhite_cvBlack:favW_tok",
                        title="Perceived Risk to Climate-Forcing Industry",
                        xlabel="Prior: pro-white favoritism in government", 
                        ylabel="Marginal effect of majority-white workforce", 
                        factor_labels=c("No","Yes")) # Figure H1, left-hand panel

lucid %>%
  filter(treatment != "control") %>%
  mutate(treatment2 = as.integer(treatment == "ffBlack_cvWhite")) %>%
  lm(decline_expectation_cv ~ treatment2*favW_tok +
       gender_fct + ideology + party_fct + age + income + education +
       as.factor(NONHWHITE),
     data = .) -> risk_intV
interaction_plot_binary(model = risk_intV,
                        effect = "treatment2",
                        moderator = "favW_tok", 
                        interaction = "treatment2:favW_tok",
                        title="Perceived Risk to Climate-Vulnerable Industry",
                        xlabel="Prior: pro-white favoritism in government", 
                        ylabel="Marginal effect of majority-white workforce", 
                        factor_labels=c("No","Yes")) # Figure H1, right-hand panel

# > Appendix I. Perceptions of Climate Risk: Excluding Inattentive Subjects ####

lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  lm_robust(decline_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg1a_attn # I1, Model 1
lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  lm_robust(decline_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg1_attn # I1, Model 2

lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(decline_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg2a_attn # I1, Model 3
lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(decline_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg2_attn # I1, Model 4

lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 0) %>%
  lm_robust(decline_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg3a_attn # I1, Model 5
lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 0) %>%
  lm_robust(decline_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg3_attn # I1, Model 6

#NOTE: the coefficients in I1 are inverse of what's shown here for Models 7-12,
#as "treatment" indicates a Black-majority climate-vulnerable industry,
#while the table lists results for the inverse (white-majority climate-vulnerable industry)
lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  lm_robust(decline_expectation_cv ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg4a_attn # I1, Model 7
lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  lm_robust(decline_expectation_cv ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg4_attn # I1, Model 8

lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(decline_expectation_cv ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg5a_attn # I1, Model 9
lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(decline_expectation_cv ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg5_attn # I1, Model 10

lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 0) %>%
  lm_robust(decline_expectation_cv ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg6a_attn # I1, Model 11
lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 0) %>%
  lm_robust(decline_expectation_cv ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_reg6_attn # I1, Model 12

# > Appendix J. Access to Subsidies: Bonferroni Correction ####

c(bias_sub_reg1$p.value[[1]],
  bias_sub_reg2$p.value[[1]],
  bias_sub_reg3$p.value[[1]]) -> sub_p_vals
p.adjust(p = sub_p_vals, method = "bonferroni")

# > Appendix K. Access to Subsidies: Randomization Inference ####

# >> All ####

lucid_ri_fav_sub <- lucid %>%
  dplyr::select(treatment, 
                subsidy_expectation_ff,
                NONHWHITE,
                gender_fct,
                ideology,
                party_fct,
                age,
                income,
                education) %>%
  mutate(treatment = case_when(treatment == "ffWhite_cvBlack" ~ 1,
                               treatment == "ffBlack_cvWhite" ~ 0,
                               treatment == "control" ~ NA_real_,
                               TRUE ~ NA_real_)) %>%
  rename(Z = treatment,
         Y = subsidy_expectation_ff,
         block = NONHWHITE) %>%
  filter(!is.na(Z) & !is.na(block) & !is.na(Y) &
           !is.na(gender_fct) & !is.na(ideology) & !is.na(party_fct) &
           !is.na(age) & !is.na(income) & !is.na(education))

with(lucid_ri_fav_sub, table(block, Z))

block_m_sub <- with(lucid_ri_fav_sub, tapply(Z, block, sum) / 2)

declaration_sub <- with(lucid_ri_fav_sub, {
  declare_ra(blocks = block,
             block_m = block_m_sub)
})

subF_ri_out <- conduct_ri(Y ~ Z +
                            gender_fct + ideology + party_fct + age + income + education,
                          sharp_hypothesis = 0,
                          declaration = declaration_sub,
                          data = lucid_ri_fav_sub[!is.na(lucid_ri_fav_sub$Z) &
                                                    !is.na(lucid_ri_fav_sub$Y), ])

tidy(subF_ri_out) # Table K1, Model 1

# >> Prior: whites favored ####

lucid_ri_fav_sub2 <- lucid %>%
  filter(favW_tok == 1) %>%
  dplyr::select(treatment, 
                subsidy_expectation_ff,
                NONHWHITE,
                gender_fct,
                ideology,
                party_fct,
                age,
                income,
                education) %>%
  mutate(treatment = case_when(treatment == "ffWhite_cvBlack" ~ 1,
                               treatment == "ffBlack_cvWhite" ~ 0,
                               treatment == "control" ~ NA_real_,
                               TRUE ~ NA_real_)) %>%
  rename(Z = treatment,
         Y = subsidy_expectation_ff,
         block = NONHWHITE) %>%
  filter(!is.na(Z) & !is.na(block) & !is.na(Y) &
           !is.na(gender_fct) & !is.na(ideology) & !is.na(party_fct) &
           !is.na(age) & !is.na(income) & !is.na(education))

with(lucid_ri_fav_sub2, table(block, Z))

block_m_sub2 <- with(lucid_ri_fav_sub2, tapply(Z, block, sum) / 2)

declaration_sub2 <- with(lucid_ri_fav_sub2, {
  declare_ra(blocks = block,
             block_m = block_m_sub2)
})

subF_ri_out2 <- conduct_ri(Y ~ Z +
                             gender_fct + ideology + party_fct + age + income + education,
                           sharp_hypothesis = 0,
                           declaration = declaration_sub2,
                           data = lucid_ri_fav_sub2[!is.na(lucid_ri_fav_sub2$Z) &
                                                      !is.na(lucid_ri_fav_sub2$Y), ])
tidy(subF_ri_out2) # Table K1, Model 2

# >> Prior: whites not favored ####

lucid_ri_fav_sub3 <- lucid %>%
  filter(favW_tok == 0) %>%
  dplyr::select(treatment, 
                subsidy_expectation_ff,
                NONHWHITE,
                gender_fct,
                ideology,
                party_fct,
                age,
                income,
                education) %>%
  mutate(treatment = case_when(treatment == "ffWhite_cvBlack" ~ 1,
                               treatment == "ffBlack_cvWhite" ~ 0,
                               treatment == "control" ~ NA_real_,
                               TRUE ~ NA_real_)) %>%
  rename(Z = treatment,
         Y = subsidy_expectation_ff,
         block = NONHWHITE) %>%
  filter(!is.na(Z) & !is.na(block) & !is.na(Y) &
           !is.na(gender_fct) & !is.na(ideology) & !is.na(party_fct) &
           !is.na(age) & !is.na(income) & !is.na(education))

with(lucid_ri_fav_sub3, table(block, Z))

block_m_sub3 <- with(lucid_ri_fav_sub3, tapply(Z, block, sum) / 2)

declaration_sub3 <- with(lucid_ri_fav_sub3, {
  declare_ra(blocks = block,
             block_m = block_m_sub3)
})

subF_ri_out3 <- conduct_ri(Y ~ Z +
                             gender_fct + ideology + party_fct + age + income + education,
                           sharp_hypothesis = 0,
                           declaration = declaration_sub3,
                           data = lucid_ri_fav_sub3[!is.na(lucid_ri_fav_sub3$Z) &
                                                      !is.na(lucid_ri_fav_sub3$Y), ])
tidy(subF_ri_out3)

# > Appendix L. Access to Subsidies: Alternative Bias Measures ####

# >> L1. White people have an easier time... ####

lucid %>%
  filter(treatment != "control") %>%
  filter(whitesupport %in% 3:4) %>%
  lm_robust(subsidy_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg2a_alt1 # Table L1, Model 1
lucid %>%
  filter(treatment != "control") %>%
  filter(whitesupport %in% 3:4) %>%
  lm_robust(subsidy_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg2_alt1 # Table L1, Model 2

lucid %>%
  filter(treatment != "control") %>%
  filter(whitesupport %in% 0:2) %>%
  lm_robust(subsidy_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg3a_alt1 # Table L1, Model 3
lucid %>%
  filter(treatment != "control") %>%
  filter(whitesupport %in% 0:2) %>%
  lm_robust(subsidy_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg3_alt1 # Table L1, Model 4

# >> L2. On average, politicians are today more responsive... ####

lucid %>%
  filter(treatment != "control") %>%
  filter(whiteresponsiveness %in% 3:4) %>%
  lm_robust(subsidy_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg2a_alt2 # Table L2, Model 1
lucid %>%
  filter(treatment != "control") %>%
  filter(whiteresponsiveness %in% 3:4) %>%
  lm_robust(subsidy_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg2_alt2 # Table L2, Model 2

lucid %>%
  filter(treatment != "control") %>%
  filter(whiteresponsiveness %in% 0:2) %>%
  lm_robust(subsidy_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg3a_alt2 # Table L2, Model 3
lucid %>%
  filter(treatment != "control") %>%
  filter(whiteresponsiveness %in% 0:2) %>%
  lm_robust(subsidy_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg3_alt2 # Table L2, Model 4

# >> L3. When the government spends money, white people... ####

lucid %>%
  filter(treatment != "control") %>%
  filter(whitemoneywins %in% 3:4) %>%
  lm_robust(subsidy_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg2a_alt3 # Table L3, Model 1
lucid %>%
  filter(treatment != "control") %>%
  filter(whitemoneywins %in% 3:4) %>%
  lm_robust(subsidy_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg2_alt3 # Table L3, Model 2

lucid %>%
  filter(treatment != "control") %>%
  filter(whitemoneywins %in% 0:2) %>%
  lm_robust(subsidy_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg3a_alt3 # Table L3, Model 3
lucid %>%
  filter(treatment != "control") %>%
  filter(whitemoneywins %in% 0:2) %>%
  lm_robust(subsidy_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg3_alt3 # Table L3, Model 4

# > Appendix M. Access to Subsidies: Interaction Model ####

lucid %>%
  filter(treatment != "control") %>%
  lm(subsidy_expectation_ff ~ treatment*favW_tok +
       gender_fct + ideology + party_fct + age + income + education +
       as.factor(NONHWHITE),
     data = .) -> sub_intF
interaction_plot_binary(model = sub_intF,
                        effect = "treatmentffWhite_cvBlack",
                        moderator = "favW_tok", 
                        interaction = "treatmentffWhite_cvBlack:favW_tok",
                        title="Pr(Subsidies for Climate Forcing Industry = 1)",
                        xlabel="Prior: pro-white favoritism in government", 
                        ylabel="Marginal effect of majority-white workforce", 
                        factor_labels=c("No","Yes")) # Figure M1

# > Appendix N. Access to Subsidies: Excluding Inattentive Subjects ####

lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  lm_robust(subsidy_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg1a_attn # N1, Model 1
lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  lm_robust(subsidy_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg1_attn # N1, Model 2

lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(subsidy_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg2a_attn # N1, Model 3
lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  lm_robust(subsidy_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg2_attn # N1, Model 4

lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 0) %>%
  lm_robust(subsidy_expectation_ff ~ treatment,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg3a_attn # N1, Model 5
lucid %>%
  filter(attncheck1 == "Climate change") %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 0) %>%
  lm_robust(subsidy_expectation_ff ~ treatment  +
              gender_fct + ideology + party_fct + age + income + education,
            data = .,
            fixed_effects = ~ NONHWHITE) -> bias_sub_reg3_attn # N1, Model 6

# > Appendix O. Mediation Analysis ####

set.seed(430)

# >> O1. Climate-forcing industry ####
lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  filter(!is.na(subsidy_expectation_ff) & 
           !is.na(treatment) &
           !is.na(decline_expectation_ff)) %>%
  lm(subsidy_expectation_ff ~ treatment + as.factor(NONHWHITE),
     data = .,
     na.action = na.omit) -> med_fav1
lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  filter(!is.na(subsidy_expectation_ff) & 
           !is.na(treatment) &
           !is.na(decline_expectation_ff)) %>%
  lm(decline_expectation_ff ~ subsidy_expectation_ff + treatment + as.factor(NONHWHITE),
     data = .,
     na.action = na.omit) -> med_fav2
mediate(med_fav1, med_fav2,
        treat = "treatment",
        mediator = "subsidy_expectation_ff") -> med_fav_out
summary(med_fav_out) # Table O1

# >> O2. Climate-vulnerable industry ####
lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  filter(!is.na(subsidy_expectation_ff) & 
           !is.na(treatment) &
           !is.na(decline_expectation_cv)) %>%
  lm(subsidy_expectation_ff ~ treatment + as.factor(NONHWHITE),
     data = .,
     na.action = na.omit) -> med_fav1a
lucid %>%
  filter(treatment != "control") %>%
  filter(favW_tok == 1) %>%
  filter(!is.na(subsidy_expectation_ff) & 
           !is.na(treatment) &
           !is.na(decline_expectation_cv)) %>%
  lm(decline_expectation_cv ~ subsidy_expectation_ff + treatment + as.factor(NONHWHITE),
     data = .,
     na.action = na.omit) -> med_fav2a
mediate(med_fav1a, med_fav2a,
        treat = "treatment",
        mediator = "subsidy_expectation_ff") -> med_fav_out2
summary(med_fav_out2) # Table O2

# > Appendix P. Analysis of Written Responses ####

# >> P2. Structural Topic Modeling ####

text <- data.frame(cbind(
  lucid$treatment,
  lucid$short_answer
))
colnames(text) <- c("treat", "text")

Corpus <- Corpus(VectorSource(text$text))
Corpus.clean <- tm_map(Corpus, stripWhitespace)
Corpus.clean <- tm_map(Corpus.clean, removePunctuation)
Corpus.clean <- tm_map(Corpus.clean, tolower)
Corpus.clean <- tm_map(Corpus.clean, removeWords, stopwords('english'))

dataframe <- data.frame(text=unlist(sapply(Corpus.clean, `[`)), 
                        stringsAsFactors=FALSE)

text <- cbind(lucid, dataframe)
text <- text %>%
  as.data.frame()

processed <- textProcessor(text$text, metadata = text)
out <- prepDocuments(processed$documents, processed$vocab, processed$meta)

set.seed(101)

text.model <- stm(documents = out$documents, vocab = out$vocab,
                  K = 4, prevalence =~ treatment,
                  max.em.its = 75, data = out$meta,
                  init.type = "Spectral")

labelTopics(text.model)
plot.STM(text.model,type = "labels", topics = c(1,2,3,4)) # Figure P1

# >> P1. Subsidies Expectations and Mentions of "White" ####

text <- text %>%
  mutate(white_in_response = str_detect(tolower(short_answer),
                                        pattern = "white"),
         black_in_response = str_detect(tolower(short_answer),
                                        pattern = "black"))

text %>%
  filter(treatment != "control") %>%
  mutate(white_sub = as.logical((subsidy_expectation_ff == 1 &  # identifying choice of subsidies for white-majority industry
                                   treatment == "ffWhite_cvBlack") | 
                                  (subsidy_expectation_ff == 0 & 
                                     treatment == "ffBlack_cvWhite"))) %>%
  feols(white_sub ~ white_in_response,
        se = "hetero",
        data = .) %>%
  tidy() # Table P1

# > Appendix Q. Worker Security and Political Efficacy: Bonferroni Correction ####

c(bias_protection_reg1$pval[[1]],
  bias_protection_reg2$pval[[1]],
  bias_protection_reg3$pval[[1]],
  bias_protection_reg4$pval[[1]],
  bias_protection_reg5$pval[[1]],
  bias_protection_reg6$pval[[1]]) -> worker_p_vals
p.adjust(p = worker_p_vals, method = "bonferroni")

# > Appendix R. Worker Security and Political Efficacy: Alternative Bias Measures ####

alt_measures_by_respondent <- lucid %>%
  dplyr::select(response_id, whitesupport, whiteresponsiveness, whitemoneywins)
lucid2 %<>% left_join(alt_measures_by_respondent, by = "response_id") # joining subject-level covariates with subject-task-level dataset

# >> R1. White people have an easier time... ####

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = whitesupport %in% 3:4,
     na.action = na.omit) -> bias_protection_reg2_alt1 # Table R1, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = whitesupport %in% 0:2,
     na.action = na.omit) -> bias_protection_reg3_alt1 # Table R1, Model 2
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = whitesupport %in% 3:4,
     na.action = na.omit) -> bias_protection_reg5_alt1 # Table R1, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = whitesupport %in% 0:2,
     na.action = na.omit) -> bias_protection_reg6_alt1 # Table R1, Model 4

# >> R2. On average, politicians are today more responsive... ####

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = whiteresponsiveness %in% 3:4,
     na.action = na.omit) -> bias_protection_reg2_alt2 # Table R2, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = whiteresponsiveness %in% 0:2,
     na.action = na.omit) -> bias_protection_reg3_alt2 # Table R2, Model 2
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = whiteresponsiveness %in% 3:4,
     na.action = na.omit) -> bias_protection_reg5_alt2 # Table R2, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = whiteresponsiveness %in% 0:2,
     na.action = na.omit) -> bias_protection_reg6_alt2 # Table R2, Model 4

# >> R3. When the government spends money, white people... ####

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = whitemoneywins %in% 3:4,
     na.action = na.omit) -> bias_protection_reg2_alt3 # Table R3, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = whitemoneywins %in% 0:2,
     na.action = na.omit) -> bias_protection_reg3_alt3 # Table R3, Model 2
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = whitemoneywins %in% 3:4,
     na.action = na.omit) -> bias_protection_reg5_alt3 # Table R3, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = whitemoneywins %in% 0:2,
     na.action = na.omit) -> bias_protection_reg6_alt3 # Table R3, Model 4

# > Appendix S. Worker Security and Political Efficacy: Excluding Inattentive Subjects ####

attncheck <- lucid[, c("response_id", "attncheck1")]
lucid2 %<>% left_join(attncheck, by = "response_id") # joining subject-level attention measure with subject-task-level dataset

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2[lucid2$attncheck1 == "Climate change", ],
     na.action = na.omit) -> bias_protection_reg1_attn # Table S1, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2[lucid2$attncheck1 == "Climate change", ],
     subset = favW_tok == 1,
     na.action = na.omit) -> bias_protection_reg2_attn # Table S1, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2[lucid2$attncheck1 == "Climate change", ],
     subset = favW_tok == 0,
     na.action = na.omit) -> bias_protection_reg3_attn # Table S1, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2[lucid2$attncheck1 == "Climate change", ],
     na.action = na.omit) -> bias_protection_reg4_attn # Table S1, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2[lucid2$attncheck1 == "Climate change", ],
     subset = favW_tok == 1,
     na.action = na.omit) -> bias_protection_reg5_attn # Table S1, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2[lucid2$attncheck1 == "Climate change", ],
     subset = favW_tok == 0,
     na.action = na.omit) -> bias_protection_reg6_attn # Table S1, Model 6

# >> First response only (paragraph in Appendix S) ####

lucid %>%
  filter(ind2x2_order_first %in% c("FFWhite", "FFBlack", "CVWhite", "CVBlack")) %>%
  feols(aid_2x2_first ~ I(ind2x2_order_first %in% c("FFWhite","CVWhite")),
        se = "hetero") %>%
  tidy()
lucid %>%
  filter(ind2x2_order_first %in% c("FFWhite", "FFBlack", "CVWhite", "CVBlack")) %>%
  feols(mob_2x2_first ~ I(ind2x2_order_first %in% c("FFWhite","CVWhite")),
        se = "hetero") %>%
  tidy()

# > Appendix T. Worker Security and Political Efficacy: By Industry ####

# >> T1. Climate-forcing industries ####

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = industry %in% c("ffblack_aid", "ffwhite_aid"),
     na.action = na.omit) -> bias_protection_reg1_ff # Table T1, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & industry %in% c("ffblack_aid", "ffwhite_aid"),
     na.action = na.omit) -> bias_protection_reg2_ff # Table T1, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & industry %in% c("ffblack_aid", "ffwhite_aid"),
     na.action = na.omit) -> bias_protection_reg3_ff # Table T1, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = industry %in% c("ffblack_mob", "ffwhite_mob"),
     na.action = na.omit) -> bias_protection_reg4_ff # Table T1, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & industry %in% c("ffblack_mob", "ffwhite_mob"),
     na.action = na.omit) -> bias_protection_reg5_ff # Table T1, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & industry %in% c("ffblack_mob", "ffwhite_mob"),
     na.action = na.omit) -> bias_protection_reg6_ff # Table T1, Model 6

# >> T2. Climate-vulnerable industries ####

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = industry %in% c("cvblack_aid", "cvwhite_aid"),
     na.action = na.omit) -> bias_protection_reg1_cv # Table T2, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & industry %in% c("cvblack_aid", "cvwhite_aid"),
     na.action = na.omit) -> bias_protection_reg2_cv # Table T2, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & industry %in% c("cvblack_aid", "cvwhite_aid"),
     na.action = na.omit) -> bias_protection_reg3_cv # Table T2, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = industry %in% c("cvblack_mob", "cvwhite_mob"),
     na.action = na.omit) -> bias_protection_reg4_cv # Table T2, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & industry %in% c("cvblack_mob", "cvwhite_mob"),
     na.action = na.omit) -> bias_protection_reg5_cv # Table T2, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & industry %in% c("cvblack_mob", "cvwhite_mob"),
     na.action = na.omit) -> bias_protection_reg6_cv # Table T2, Model 6

# > Appendix U. Worker Security and Political Efficacy: Other Subgroups ####

# >> U1. Democrats vs. Republicans ####

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = party_bin == "Democrat",
     na.action = na.omit) -> bias_protection_reg1_dem # Table U1, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & party_bin == "Democrat",
     na.action = na.omit) -> bias_protection_reg2_dem # Table U1, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & party_bin == "Democrat",
     na.action = na.omit) -> bias_protection_reg3_dem # Table U1, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = party_bin == "Democrat",
     na.action = na.omit) -> bias_protection_reg4_dem # Table U1, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & party_bin == "Democrat",
     na.action = na.omit) -> bias_protection_reg5_dem # Table U1, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & party_bin == "Democrat",
     na.action = na.omit) -> bias_protection_reg6_dem # Table U1, Model 6

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = party_bin == "Republican",
     na.action = na.omit) -> bias_protection_reg1_gop # Tablel U2, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & party_bin == "Republican",
     na.action = na.omit) -> bias_protection_reg2_gop # Tablel U2, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & party_bin == "Republican",
     na.action = na.omit) -> bias_protection_reg3_gop # Tablel U2, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = party_bin == "Republican",
     na.action = na.omit) -> bias_protection_reg4_gop # Tablel U2, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & party_bin == "Republican",
     na.action = na.omit) -> bias_protection_reg5_gop # Tablel U2, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & party_bin == "Republican",
     na.action = na.omit) -> bias_protection_reg6_gop # Tablel U2, Model 6

# >> U2. Biden Voters vs. Trump Voters ####

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = presvote == "Joe Biden",
     na.action = na.omit) -> bias_protection_reg1_biden # Table U3, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & presvote == "Joe Biden",
     na.action = na.omit) -> bias_protection_reg2_biden # Table U3, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & presvote == "Joe Biden",
     na.action = na.omit) -> bias_protection_reg3_biden # Table U3, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = presvote == "Joe Biden",
     na.action = na.omit) -> bias_protection_reg4_biden # Table U3, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & presvote == "Joe Biden",
     na.action = na.omit) -> bias_protection_reg5_biden # Table U3, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & presvote == "Joe Biden",
     na.action = na.omit) -> bias_protection_reg6_biden # Table U3, Model 6

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = presvote == "Donald Trump",
     na.action = na.omit) -> bias_protection_reg1_trump # Table U4, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & presvote == "Donald Trump",
     na.action = na.omit) -> bias_protection_reg2_trump # Table U4, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & presvote == "Donald Trump",
     na.action = na.omit) -> bias_protection_reg3_trump # Table U4, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = presvote == "Donald Trump",
     na.action = na.omit) -> bias_protection_reg4_trump # Table U4, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & presvote == "Donald Trump",
     na.action = na.omit) -> bias_protection_reg5_trump # Table U4, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & presvote == "Donald Trump",
     na.action = na.omit) -> bias_protection_reg6_trump # Table U4, Model 6

# >> U3. Wealthy vs. Poor ####

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = income >= 3,
     na.action = na.omit) -> bias_protection_reg1_wealthy # Table U5, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & income >= 3,
     na.action = na.omit) -> bias_protection_reg2_wealthy # Table U5, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & income >= 3,
     na.action = na.omit) -> bias_protection_reg3_wealthy # Table U5, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = income >= 3,
     na.action = na.omit) -> bias_protection_reg4_wealthy # Table U5, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & income >= 3,
     na.action = na.omit) -> bias_protection_reg5_wealthy # Table U5, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & income >= 3,
     na.action = na.omit) -> bias_protection_reg6_wealthy # Table U5, Model 6

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = income == 0,
     na.action = na.omit) -> bias_protection_reg1_poor # Table U6, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & income == 0,
     na.action = na.omit) -> bias_protection_reg2_poor # Table U6, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & income == 0,
     na.action = na.omit) -> bias_protection_reg3_poor # Table U6, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = income == 0,
     na.action = na.omit) -> bias_protection_reg4_poor # Table U6, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & income == 0,
     na.action = na.omit) -> bias_protection_reg5_poor # Table U6, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & income == 0,
     na.action = na.omit) -> bias_protection_reg6_poor # Table U6, Model 6

# >> U4. Racial Resentment ####

alt_measures_by_respondent2 <- lucid %>%
  dplyr::select(response_id, INDEX_animus_reverse)
lucid2 %<>% left_join(alt_measures_by_respondent2, by = "response_id")

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = INDEX_animus_reverse >= median(lucid2$INDEX_animus_reverse, na.rm = TRUE),
     na.action = na.omit) -> bias_protection_reg1_resent # Table U7, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & INDEX_animus_reverse >= median(lucid2$INDEX_animus_reverse, na.rm = TRUE),
     na.action = na.omit) -> bias_protection_reg2_resent # Table U7, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & INDEX_animus_reverse >= median(lucid2$INDEX_animus_reverse, na.rm = TRUE),
     na.action = na.omit) -> bias_protection_reg3_resent # Table U7, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = INDEX_animus_reverse >= median(lucid2$INDEX_animus_reverse, na.rm = TRUE),
     na.action = na.omit) -> bias_protection_reg4_resent # Table U7, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & INDEX_animus_reverse >= median(lucid2$INDEX_animus_reverse, na.rm = TRUE),
     na.action = na.omit) -> bias_protection_reg5_resent # Table U7, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & INDEX_animus_reverse >= median(lucid2$INDEX_animus_reverse, na.rm = TRUE),
     na.action = na.omit) -> bias_protection_reg6_resent # Table U7, Model 6

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = INDEX_animus_reverse < median(lucid2$INDEX_animus_reverse, na.rm = TRUE),
     na.action = na.omit) -> bias_protection_reg1_resentNo # Table U8, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & INDEX_animus_reverse < median(lucid2$INDEX_animus_reverse, na.rm = TRUE),
     na.action = na.omit) -> bias_protection_reg2_resentNo # Table U8, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & INDEX_animus_reverse < median(lucid2$INDEX_animus_reverse, na.rm = TRUE),
     na.action = na.omit) -> bias_protection_reg3_resentNo # Table U8, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = INDEX_animus_reverse < median(lucid2$INDEX_animus_reverse, na.rm = TRUE),
     na.action = na.omit) -> bias_protection_reg4_resentNo # Table U8, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & INDEX_animus_reverse < median(lucid2$INDEX_animus_reverse, na.rm = TRUE),
     na.action = na.omit) -> bias_protection_reg5_resentNo # Table U8, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & INDEX_animus_reverse < median(lucid2$INDEX_animus_reverse, na.rm = TRUE),
     na.action = na.omit) -> bias_protection_reg6_resentNo # Table U8, Model 6

# >> U5. Race and Ethnicity ####

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = NONHWHITE == 1,
     na.action = na.omit) -> bias_protection_reg1_nonhispwhite # Table U9, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & NONHWHITE == 1,
     na.action = na.omit) -> bias_protection_reg2_nonhispwhite # Table U9, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & NONHWHITE == 1,
     na.action = na.omit) -> bias_protection_reg3_nonhispwhite # Table U9, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = NONHWHITE == 1,
     na.action = na.omit) -> bias_protection_reg4_nonhispwhite # Table U9, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & NONHWHITE == 1,
     na.action = na.omit) -> bias_protection_reg5_nonhispwhite # Table U9, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & NONHWHITE == 1,
     na.action = na.omit) -> bias_protection_reg6_nonhispwhite # Table U9, Model 6

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = NONHWHITE == 0,
     na.action = na.omit) -> bias_protection_reg1_minority # Table U10, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & NONHWHITE == 0,
     na.action = na.omit) -> bias_protection_reg2_minority # Table U10, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & NONHWHITE == 0,
     na.action = na.omit) -> bias_protection_reg3_minority # Table U10, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = NONHWHITE == 0,
     na.action = na.omit) -> bias_protection_reg4_minority # Table U10, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & NONHWHITE == 0,
     na.action = na.omit) -> bias_protection_reg5_minority # Table U10, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & NONHWHITE == 0,
     na.action = na.omit) -> bias_protection_reg6_minority # Table U10, Model 6

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = race == "Black or African American",
     na.action = na.omit) -> bias_protection_reg1_black # Table U11, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & race == "Black or African American",
     na.action = na.omit) -> bias_protection_reg2_black # Table U11, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & race == "Black or African American",
     na.action = na.omit) -> bias_protection_reg3_black # Table U11, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = race == "Black or African American",
     na.action = na.omit) -> bias_protection_reg4_black # Table U11, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & race == "Black or African American",
     na.action = na.omit) -> bias_protection_reg5_black # Table U11, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & race == "Black or African American",
     na.action = na.omit) -> bias_protection_reg6_black # Table U11, Model 6

# >> U6. Climate Concerns ####

climateconcern <- lucid[, c("response_id", "climateconcern")]
lucid2 %<>% left_join(climateconcern, by = "response_id")

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = climateconcern %in% c("Somewhat worried", "Very worried"),
     na.action = na.omit) -> bias_protection_reg1_concerned # Table U12, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & climateconcern %in% c("Somewhat worried", "Very worried"),
     na.action = na.omit) -> bias_protection_reg2_concerned # Table U12, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & climateconcern %in% c("Somewhat worried", "Very worried"),
     na.action = na.omit) -> bias_protection_reg3_concerned # Table U12, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = climateconcern %in% c("Somewhat worried", "Very worried"),
     na.action = na.omit) -> bias_protection_reg4_concerned # Table U12, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & climateconcern %in% c("Somewhat worried", "Very worried"),
     na.action = na.omit) -> bias_protection_reg5_concerned # Table U12, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & climateconcern %in% c("Somewhat worried", "Very worried"),
     na.action = na.omit) -> bias_protection_reg6_concerned # Table U12, Model 6

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = climateconcern %in% c("Not at all worried", "Not very worried"),
     na.action = na.omit) -> bias_protection_reg1_concernedNot # Table U13, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & climateconcern %in% c("Not at all worried", "Not very worried"),
     na.action = na.omit) -> bias_protection_reg2_concernedNot # Table U13, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & climateconcern %in% c("Not at all worried", "Not very worried"),
     na.action = na.omit) -> bias_protection_reg3_concernedNot # Table U13, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = climateconcern %in% c("Not at all worried", "Not very worried"),
     na.action = na.omit) -> bias_protection_reg4_concernedNot # Table U13, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & climateconcern %in% c("Not at all worried", "Not very worried"),
     na.action = na.omit) -> bias_protection_reg5_concernedNot # Table U13, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & climateconcern %in% c("Not at all worried", "Not very worried"),
     na.action = na.omit) -> bias_protection_reg6_concernedNot # Table U13, Model 6

# >> U7. Status Threat ####

statusthreat <- lucid[, c("response_id", "raceimportant", "wayoflife", "discrim_white")]
lucid2 %<>% left_join(statusthreat, by = "response_id")

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = wayoflife %in% c("Somewhat worried", "Very worried"),
     na.action = na.omit) -> bias_protection_reg1_threatYes # Table U14, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & wayoflife %in% c("Somewhat worried", "Very worried"),
     na.action = na.omit) -> bias_protection_reg2_threatYes # Table U14, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & wayoflife %in% c("Somewhat worried", "Very worried"),
     na.action = na.omit) -> bias_protection_reg3_threatYes # Table U14, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = wayoflife %in% c("Somewhat worried", "Very worried"),
     na.action = na.omit) -> bias_protection_reg4_threatYes # Table U14, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & wayoflife %in% c("Somewhat worried", "Very worried"),
     na.action = na.omit) -> bias_protection_reg5_threatYes # Table U14, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & wayoflife %in% c("Somewhat worried", "Very worried"),
     na.action = na.omit) -> bias_protection_reg6_threatYes # Table U14, Model 6

felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = wayoflife %in% c("Not at all worried", "Not very worried"),
     na.action = na.omit) -> bias_protection_reg1_threatNo # Table U15, Model 1
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & wayoflife %in% c("Not at all worried", "Not very worried"),
     na.action = na.omit) -> bias_protection_reg2_threatNo # Table U15, Model 2
felm(answer ~ race_bin_aid | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & wayoflife %in% c("Not at all worried", "Not very worried"),
     na.action = na.omit) -> bias_protection_reg3_threatNo # Table U15, Model 3
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = wayoflife %in% c("Not at all worried", "Not very worried"),
     na.action = na.omit) -> bias_protection_reg4_threatNo # Table U15, Model 4
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 1 & wayoflife %in% c("Not at all worried", "Not very worried"),
     na.action = na.omit) -> bias_protection_reg5_threatNo # Table U15, Model 5
felm(answer ~ race_bin_mob | response_id | 0 | response_id,
     data = lucid2,
     subset = favW_tok == 0 & wayoflife %in% c("Not at all worried", "Not very worried"),
     na.action = na.omit) -> bias_protection_reg6_threatNo # Table U15, Model 6

